library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.1
library(lavaan)
## Warning: package 'lavaan' was built under R version 4.1.1
library(DiagrammeR)
## Warning: package 'DiagrammeR' was built under R version 4.1.1
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.1
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.1
# For saving SEM diagrams:
library(purrr)
## Warning: package 'purrr' was built under R version 4.1.1
library(DiagrammeRsvg)
## Warning: package 'DiagrammeRsvg' was built under R version 4.1.3
library(rsvg)
## Warning: package 'rsvg' was built under R version 4.1.3
library(png)
## Warning: package 'png' was built under R version 4.1.1
library(grid)
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.1.3

Import data

combined=read.csv("data/annual_averages/annual_data_compiled_regions.csv",stringsAsFactors = F)
cnames=read.csv("analysis/column_names_region.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975) %>% arrange(Region,Year)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
regions=unique(focaldata$region)
regionorder=c("West","North","South")
years=1975:2021

focaldata = focaldata %>% 
  mutate(tzoop=hcope+clad+mysid+pcope+rotif_m,
         tzoop_e=hcope_e+clad_e+mysid_e+pcope_e+rotif_e,
         hzoop=hcope+clad+rotif_m,
         hzoop_e=hcope_e+clad_e+rotif_e,
         pzoop=mysid+pcope,
         pzoop_e=mysid_e+pcope_e) 
fvars=c(fvars,"tzoop","tzoop_e",
        "hzoop","hzoop_e",
        "pzoop","pzoop_e")
cnames=rbind(cnames,data.frame(Longname=NA,Shortname=c("tzoop","tzoop_e",
                                                       "hzoop","hzoop_e",
                                                       "pzoop","pzoop_e"),
                               Diagramname=c("Total Zooplankton\nBiomass",
                                             "Total Zooplankton\nEnergy",
                                             "Herbivorous Zooplankton\nBiomass",
                                             "Herbivorous Zooplankton\nEnergy",
                                             "Predatory Zooplankton\nBiomass",
                                             "Predatory Zooplankton\nEnergy"),
                               Datacolumn=NA,Log="yes",
                               Color=c("black","black","#ED7D31","#ED7D31","#7030A0",
                                       "#7030A0")))

#focal variables
varnames=c("temp","flow","secchi","chla","hzoop","pzoop","tzoop","amphi","potam","corbic","estfish","estfish_bsmt","estfish_stn")

source("analysis/myLavaanPlot.r")
source("analysis/semDiagramFunctions.r")

Data prep

Log transform, scale

#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
  x2=x[which(!is.na(x))]
  if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
  else {log(x)}
}
focaldatalog = focaldata %>% 
  mutate_at(logvars,logtrans)

#scale data
fdr0=focaldatalog
tvars=fvars[-(1:2)]

fdr=fdr0 %>% group_by(region) %>% 
  #lag
  #mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

#detrended data
fdr_dtr=fdr0 %>% group_by(region) %>% 
  #detrend
  mutate_at(tvars,function(x) { 
    x<<-x
    if(!all(is.na(x))) {
      if((length(which(x==0))/length(x))<0.5) {
        x2<<-x
        x2[x2==0]=NA
        res<<-residuals(lm(x2~years))
        out=x
        out[which(!is.na(x2))]=res
        return(out)
      } else {return(x)}
    } else {return(x)}
  }) %>%
  #lag
  #mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

Time series plots

## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(varnames)` instead of `varnames` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

Other useful plots

Breakdown of total zooplankton biomass.

## Warning: Removed 15 rows containing missing values (position_stack).

Similarity of fish indices.

## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 7 row(s) containing missing values (geom_path).

Correlation between biomass and energy.

for(i in 1:length(regions)) {
  dtemp=filter(fdr,region==regions[i])
  print(regions[i])
  print(cor(dtemp$tzoop,dtemp$tzoop_e,use = "p"))
  print(cor(dtemp$hzoop,dtemp$hzoop_e,use = "p"))
  print(cor(dtemp$pzoop,dtemp$pzoop_e,use = "p"))
}
## [1] "North"
##           [,1]
## [1,] 0.9989681
##           [,1]
## [1,] 0.9958644
##           [,1]
## [1,] 0.9990414
## [1] "South"
##           [,1]
## [1,] 0.9987488
##           [,1]
## [1,] 0.9977845
##           [,1]
## [1,] 0.9993711
## [1] "West"
##          [,1]
## [1,] 0.999451
##           [,1]
## [1,] 0.9986605
##           [,1]
## [1,] 0.9994171

Cross-correlation matrices

(only sig correlations shown… no correction for multiple comparisons)

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

## Warning: `expand_scale()` is deprecated; use `expansion()` instead.

Note correlation of fish indices, or lack thereof.

SEM model

With and without detrending.

West

#1
# modwest='zoop=~hcope+mysid
#         fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         fish~zoop+flow
# '
#2
# modwest='chla~potam+flow
#         tzoop~chla+potam+flow
#         estfish_bsmt~tzoop+flow
#         estfish_bsot~tzoop+flow
# '
#3
# modwest='chla~potam+flow+temp+secchi
#         tzoop~chla+potam+flow+temp+secchi
#         amphi~chla+potam+flow+temp+secchi
#         estfish_bsmt~tzoop+amphi+flow
#         #estfish_bsot~tzoop+amphi+flow+temp+secchi
#         amphi~~tzoop
# '
#4
# modwest='chla~potam+flow+temp+secchi
#         hzoop~chla+potam+flow+temp+secchi
#         pzoop~chla+potam+flow+temp+secchi+hzoop
#         amphi~chla+potam+flow+temp+secchi
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+secchi
#         amphi~~hzoop+pzoop
# '
#5
modwest='chla~potam+flow+temp+secchi
        hzoop~chla+potam+flow+temp+secchi
        pzoop~chla+potam+flow+temp+secchi+hzoop
        fish~hzoop+pzoop+potam+flow+temp+secchi
        fish=~estfish+estfish_stn+estfish_bsmt
'

modfitwest=sem(modwest, data=filter(fdr,region=="West"))
modfitwest_dtr=sem(modwest, data=filter(fdr_dtr,region=="West"))
summary(modfitwest, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 37 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                18.963
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.215
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.881    0.885
##     estfish_stn       0.907    0.113    8.042    0.000    0.799    0.877
##     estfish_bsmt      0.844    0.140    6.018    0.000    0.744    0.753
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam            -0.303    0.171   -1.775    0.076   -0.303   -0.303
##     flow              0.187    0.171    1.095    0.273    0.187    0.193
##     temp             -0.290    0.152   -1.905    0.057   -0.290   -0.287
##     secchi            0.169    0.198    0.851    0.395    0.169    0.163
##   hzoop ~                                                               
##     chla              0.511    0.097    5.260    0.000    0.511    0.589
##     potam            -0.267    0.109   -2.455    0.014   -0.267   -0.308
##     flow              0.233    0.107    2.186    0.029    0.233    0.277
##     temp              0.018    0.098    0.183    0.855    0.018    0.020
##     secchi            0.278    0.123    2.265    0.023    0.278    0.309
##   pzoop ~                                                               
##     chla              0.342    0.088    3.895    0.000    0.342    0.393
##     potam            -0.013    0.081   -0.158    0.874   -0.013   -0.015
##     flow             -0.291    0.078   -3.717    0.000   -0.291   -0.345
##     temp             -0.002    0.068   -0.030    0.976   -0.002   -0.002
##     secchi           -0.088    0.091   -0.971    0.332   -0.088   -0.097
##     hzoop             0.651    0.110    5.920    0.000    0.651    0.648
##   fish ~                                                                
##     hzoop             0.152    0.161    0.945    0.345    0.172    0.149
##     pzoop             0.343    0.145    2.366    0.018    0.390    0.340
##     potam            -0.317    0.089   -3.562    0.000   -0.359   -0.359
##     flow              0.225    0.096    2.348    0.019    0.256    0.264
##     temp              0.005    0.071    0.072    0.942    0.006    0.006
##     secchi           -0.199    0.099   -2.016    0.044   -0.226   -0.218
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.215    0.062    3.486    0.000    0.215    0.217
##    .estfish_stn       0.191    0.053    3.574    0.000    0.191    0.230
##    .estfish_bsmt      0.422    0.101    4.162    0.000    0.422    0.433
##    .chla              0.738    0.165    4.472    0.000    0.738    0.739
##    .hzoop             0.279    0.062    4.472    0.000    0.279    0.370
##    .pzoop             0.135    0.030    4.472    0.000    0.135    0.178
##    .fish              0.059    0.039    1.519    0.129    0.076    0.076
## 
## R-Square:
##                    Estimate
##     estfish           0.783
##     estfish_stn       0.770
##     estfish_bsmt      0.567
##     chla              0.261
##     hzoop             0.630
##     pzoop             0.822
##     fish              0.924
summary(modfitwest_dtr, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 32 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                17.630
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.283
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.679    0.736
##     estfish_stn       1.118    0.228    4.896    0.000    0.759    0.771
##     estfish_bsmt      0.924    0.232    3.986    0.000    0.628    0.636
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam             0.067    0.161    0.419    0.675    0.067    0.068
##     flow              0.357    0.181    1.973    0.049    0.357    0.367
##     temp             -0.350    0.154   -2.274    0.023   -0.350   -0.344
##     secchi            0.270    0.187    1.444    0.149    0.270    0.267
##   hzoop ~                                                               
##     chla              0.509    0.120    4.223    0.000    0.509    0.509
##     potam            -0.054    0.123   -0.436    0.663   -0.054   -0.054
##     flow              0.437    0.145    3.019    0.003    0.437    0.449
##     temp             -0.023    0.125   -0.185    0.854   -0.023   -0.023
##     secchi            0.439    0.146    3.004    0.003    0.439    0.435
##   pzoop ~                                                               
##     chla              0.389    0.095    4.113    0.000    0.389    0.461
##     potam             0.055    0.080    0.689    0.491    0.055    0.067
##     flow             -0.196    0.105   -1.878    0.060   -0.196   -0.240
##     temp             -0.039    0.081   -0.474    0.635   -0.039   -0.045
##     secchi            0.047    0.106    0.441    0.660    0.047    0.055
##     hzoop             0.422    0.103    4.092    0.000    0.422    0.501
##   fish ~                                                                
##     hzoop             0.038    0.117    0.330    0.742    0.057    0.056
##     pzoop             0.282    0.131    2.153    0.031    0.415    0.349
##     potam            -0.157    0.079   -1.978    0.048   -0.231   -0.234
##     flow              0.444    0.116    3.823    0.000    0.654    0.671
##     temp             -0.048    0.077   -0.625    0.532   -0.070   -0.069
##     secchi           -0.044    0.101   -0.442    0.658   -0.066   -0.065
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.390    0.101    3.870    0.000    0.390    0.458
##    .estfish_stn       0.394    0.107    3.670    0.000    0.394    0.406
##    .estfish_bsmt      0.581    0.139    4.178    0.000    0.581    0.596
##    .chla              0.802    0.179    4.472    0.000    0.802    0.807
##    .hzoop             0.466    0.104    4.472    0.000    0.466    0.469
##    .pzoop             0.198    0.044    4.472    0.000    0.198    0.281
##    .fish              0.041    0.048    0.867    0.386    0.090    0.090
## 
## R-Square:
##                    Estimate
##     estfish           0.542
##     estfish_stn       0.594
##     estfish_bsmt      0.404
##     chla              0.193
##     hzoop             0.531
##     pzoop             0.719
##     fish              0.910
labelswest <- createLabels(modfitwest, cnames)

# residuals(modfitwest)
# modificationindices(modfitwest)

North

#1
# modnorth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '
# modnorth='zoop=~clad
#         zoop~chla+corbic+potam+flow
#         chla~corbic+potam+flow
#         estfish_bsmt~zoop+flow+sside+chla
#         estfish_bsot~zoop+flow+sside+chla
# '
#2
# modnorth='chla~corbic+potam+flow
#         tzoop~chla+corbic+potam+flow
#         estfish_bsmt~tzoop+flow+chla
#         estfish_bsot~tzoop+flow
# '
#3
# modnorth='chla~corbic+potam+flow+temp+secchi
#         tzoop~chla+corbic+potam+flow+temp+secchi
#         amphi~chla+corbic+potam+flow+temp+secchi
#         estfish_bsmt~tzoop+amphi+flow+temp+secchi+chla+sside
#         #estfish_bsot~tzoop+amphi+flow+temp+secchi+sside
#         amphi~~tzoop
# '
#4
# modnorth='chla~corbic+flow+temp+secchi
#         hzoop~chla+corbic+flow+temp+secchi
#         pzoop~chla+corbic+flow+temp+secchi+hzoop
#         amphi~chla+corbic+flow+temp+secchi
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+secchi+chla+sside
#         amphi~~hzoop+pzoop
# '
#5
modnorth='chla~corbic+flow+temp+secchi
        hzoop~chla+corbic+flow+temp+secchi
        pzoop~chla+corbic+flow+temp+secchi+hzoop
        fish~chla+hzoop+pzoop+corbic+flow+temp+secchi
        fish=~estfish+estfish_stn #+estfish_bsmt
        estfish_stn~~chla
'

modfitnorth=sem(modnorth, data=filter(fdr,region=="North"))
modfitnorth_dtr=sem(modnorth, data=filter(fdr_dtr,region=="North"))
summary(modfitnorth, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 35 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 5.374
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.372
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.829    0.809
##     estfish_stn       0.900    0.175    5.153    0.000    0.747    0.844
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.465    0.129    3.589    0.000    0.465    0.466
##     flow              0.128    0.139    0.917    0.359    0.128    0.128
##     temp              0.153    0.131    1.172    0.241    0.153    0.149
##     secchi           -0.129    0.135   -0.960    0.337   -0.129   -0.126
##   hzoop ~                                                               
##     chla              0.125    0.119    1.053    0.292    0.125    0.145
##     corbic            0.313    0.122    2.569    0.010    0.313    0.364
##     flow             -0.028    0.118   -0.237    0.813   -0.028   -0.032
##     temp              0.042    0.111    0.374    0.708    0.042    0.047
##     secchi           -0.396    0.108   -3.676    0.000   -0.396   -0.448
##   pzoop ~                                                               
##     chla              0.534    0.087    6.115    0.000    0.534    0.557
##     corbic            0.373    0.095    3.915    0.000    0.373    0.390
##     flow             -0.466    0.086   -5.449    0.000   -0.466   -0.486
##     temp              0.028    0.081    0.346    0.729    0.028    0.028
##     secchi           -0.110    0.090   -1.221    0.222   -0.110   -0.112
##     hzoop             0.178    0.113    1.568    0.117    0.178    0.160
##   fish ~                                                                
##     chla             -0.505    0.181   -2.792    0.005   -0.609   -0.617
##     hzoop             0.456    0.159    2.872    0.004    0.550    0.481
##     pzoop             0.481    0.208    2.311    0.021    0.580    0.564
##     corbic           -0.248    0.157   -1.579    0.114   -0.299   -0.304
##     flow              0.147    0.150    0.981    0.327    0.177    0.180
##     temp              0.077    0.110    0.698    0.485    0.092    0.091
##     secchi           -0.306    0.125   -2.455    0.014   -0.369   -0.364
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_stn ~~                                                        
##    .chla              0.280    0.115    2.439    0.015    0.280    0.586
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.364    0.117    3.112    0.002    0.364    0.346
##    .estfish_stn       0.317    0.109    2.905    0.004    0.317    0.404
##    .chla              0.722    0.159    4.528    0.000    0.722    0.702
##    .hzoop             0.418    0.092    4.528    0.000    0.418    0.546
##    .pzoop             0.220    0.049    4.528    0.000    0.220    0.233
##    .fish              0.213    0.097    2.195    0.028    0.310    0.310
## 
## R-Square:
##                    Estimate
##     estfish           0.654
##     estfish_stn       0.596
##     chla              0.298
##     hzoop             0.454
##     pzoop             0.767
##     fish              0.690
summary(modfitnorth_dtr, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 37 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 6.333
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.275
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.595    0.600
##     estfish_stn       1.228    0.476    2.582    0.010    0.731    0.732
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.545    0.157    3.465    0.001    0.545    0.514
##     flow              0.071    0.148    0.484    0.628    0.071    0.071
##     temp              0.146    0.132    1.106    0.269    0.146    0.142
##     secchi           -0.234    0.139   -1.687    0.092   -0.234   -0.227
##   hzoop ~                                                               
##     chla              0.354    0.165    2.139    0.032    0.354    0.360
##     corbic           -0.031    0.197   -0.156    0.876   -0.031   -0.029
##     flow              0.141    0.173    0.814    0.415    0.141    0.142
##     temp              0.083    0.156    0.534    0.593    0.083    0.082
##     secchi            0.001    0.167    0.008    0.993    0.001    0.001
##   pzoop ~                                                               
##     chla              0.704    0.105    6.690    0.000    0.704    0.688
##     corbic            0.307    0.119    2.585    0.010    0.307    0.283
##     flow             -0.514    0.105   -4.884    0.000   -0.514   -0.497
##     temp              0.046    0.094    0.484    0.628    0.046    0.043
##     secchi            0.016    0.101    0.156    0.876    0.016    0.015
##     hzoop             0.040    0.094    0.420    0.675    0.040    0.038
##   fish ~                                                                
##     chla             -0.336    0.212   -1.584    0.113   -0.565   -0.576
##     hzoop             0.078    0.102    0.759    0.448    0.131    0.131
##     pzoop             0.259    0.181    1.430    0.153    0.435    0.454
##     corbic           -0.354    0.192   -1.848    0.065   -0.595   -0.573
##     flow              0.170    0.157    1.083    0.279    0.285    0.288
##     temp              0.098    0.114    0.863    0.388    0.165    0.163
##     secchi           -0.092    0.123   -0.750    0.453   -0.154   -0.153
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_stn ~~                                                        
##    .chla              0.457    0.188    2.432    0.015    0.457    0.682
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.630    0.177    3.559    0.000    0.630    0.640
##    .estfish_stn       0.600    0.249    2.406    0.016    0.600    0.602
##    .chla              0.749    0.165    4.528    0.000    0.749    0.719
##    .hzoop             0.838    0.185    4.528    0.000    0.838    0.834
##    .pzoop             0.306    0.068    4.528    0.000    0.306    0.281
##    .fish              0.182    0.117    1.560    0.119    0.513    0.513
## 
## R-Square:
##                    Estimate
##     estfish           0.360
##     estfish_stn       0.398
##     chla              0.281
##     hzoop             0.166
##     pzoop             0.719
##     fish              0.487
labelsnorth <- createLabels(modfitnorth, cnames)

# residuals(modfitnorth)
# modificationindices(modfitnorth)

South

#1
# modsouth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+corbic+flow
#         chla~corbic+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '
#2
# modsouth='chla~corbic+flow
#         tzoop~chla+corbic+flow
#         estfish_bsmt~tzoop+flow+corbic+sside
#         estfish_bsot~tzoop+flow+corbic+sside
# '
#3
# modsouth='chla~corbic+flow+temp+secchi
#         tzoop~chla+corbic+flow+temp+secchi
#         amphi~chla+corbic+flow+temp+secchi
#         estfish_bsmt~tzoop+amphi+flow+temp+secchi+corbic+sside
#         #estfish_bsot~tzoop+amphi+flow+temp+secchi+corbic+sside
#         amphi~~tzoop
# '
#4
# modsouth='chla~corbic+flow+temp+secchi
#         hzoop~chla+corbic+flow+temp+secchi
#         pzoop~chla+corbic+flow+temp+secchi+hzoop
#         amphi~chla+corbic+flow+temp+secchi
#         estfish_bsmt~hzoop+pzoop+amphi+flow+temp+secchi+corbic+sside
#         amphi~~hzoop+pzoop
# '
#5
modsouth='chla~corbic+flow+temp+secchi
        hzoop~chla+corbic+flow+temp+secchi
        pzoop~chla+corbic+flow+temp+secchi+hzoop
        fish~hzoop+pzoop+corbic+flow+temp+secchi
        fish=~estfish+estfish_stn+estfish_bsmt
'
modsouth_dtr='chla~corbic+flow+temp+secchi
        hzoop~chla+corbic+flow+temp+secchi
        pzoop~chla+corbic+flow+temp+secchi+hzoop
        fish~hzoop+pzoop+corbic+flow+temp+secchi
        fish=~estfish #+estfish_stn+estfish_bsmt
'
modfitsouth=sem(modsouth, data=filter(fdr,region=="South"))
modfitsouth_dtr=sem(modsouth_dtr, data=filter(fdr_dtr,region=="South"))
summary(modfitsouth, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        30
##                                                       
##                                                   Used       Total
##   Number of observations                            40          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                17.267
##   Degrees of freedom                                15
##   P-value (Chi-square)                           0.303
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.911    0.919
##     estfish_stn       0.776    0.125    6.221    0.000    0.707    0.762
##     estfish_bsmt      0.639    0.152    4.203    0.000    0.582    0.590
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.279    0.153    1.829    0.067    0.279    0.278
##     flow             -0.042    0.145   -0.287    0.774   -0.042   -0.043
##     temp              0.093    0.144    0.646    0.518    0.093    0.093
##     secchi           -0.380    0.161   -2.363    0.018   -0.380   -0.361
##   hzoop ~                                                               
##     chla              0.452    0.110    4.113    0.000    0.452    0.512
##     corbic            0.341    0.110    3.082    0.002    0.341    0.384
##     flow             -0.084    0.101   -0.832    0.405   -0.084   -0.099
##     temp             -0.075    0.101   -0.743    0.457   -0.075   -0.085
##     secchi            0.012    0.119    0.104    0.917    0.012    0.013
##   pzoop ~                                                               
##     chla              0.526    0.125    4.204    0.000    0.526    0.550
##     corbic           -0.008    0.117   -0.065    0.948   -0.008   -0.008
##     flow             -0.101    0.097   -1.043    0.297   -0.101   -0.110
##     temp              0.043    0.096    0.450    0.653    0.043    0.045
##     secchi           -0.185    0.114   -1.625    0.104   -0.185   -0.184
##     hzoop             0.253    0.151    1.679    0.093    0.253    0.234
##   fish ~                                                                
##     hzoop             0.220    0.125    1.762    0.078    0.242    0.213
##     pzoop            -0.090    0.110   -0.819    0.413   -0.099   -0.095
##     corbic            0.189    0.098    1.922    0.055    0.207    0.207
##     flow             -0.096    0.082   -1.173    0.241   -0.105   -0.110
##     temp              0.013    0.081    0.163    0.870    0.014    0.014
##     secchi           -0.767    0.102   -7.530    0.000   -0.842   -0.800
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.153    0.069    2.219    0.027    0.153    0.156
##    .estfish_stn       0.362    0.090    4.000    0.000    0.362    0.420
##    .estfish_bsmt      0.636    0.148    4.294    0.000    0.636    0.652
##    .chla              0.737    0.165    4.472    0.000    0.737    0.738
##    .hzoop             0.356    0.080    4.472    0.000    0.356    0.457
##    .pzoop             0.324    0.072    4.472    0.000    0.324    0.354
##    .fish              0.114    0.065    1.749    0.080    0.137    0.137
## 
## R-Square:
##                    Estimate
##     estfish           0.844
##     estfish_stn       0.580
##     estfish_bsmt      0.348
##     chla              0.262
##     hzoop             0.543
##     pzoop             0.646
##     fish              0.863
summary(modfitsouth_dtr, standardized=T, rsq=T)
## lavaan 0.6-9 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        25
##                                                       
##                                                   Used       Total
##   Number of observations                            41          47
##                                                                   
## Model Test User Model:
##                                                       
##   Test statistic                                 0.284
##   Degrees of freedom                                 1
##   P-value (Chi-square)                           0.594
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   fish =~                                                               
##     estfish           1.000                               0.909    1.000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.095    0.160    0.590    0.555    0.095    0.095
##     flow              0.063    0.166    0.378    0.706    0.063    0.065
##     temp              0.203    0.161    1.259    0.208    0.203    0.203
##     secchi            0.079    0.163    0.486    0.627    0.079    0.079
##   hzoop ~                                                               
##     chla              0.446    0.124    3.584    0.000    0.446    0.459
##     corbic            0.308    0.128    2.405    0.016    0.308    0.317
##     flow             -0.077    0.133   -0.582    0.561   -0.077   -0.083
##     temp             -0.048    0.131   -0.366    0.715   -0.048   -0.049
##     secchi            0.112    0.130    0.862    0.389    0.112    0.116
##   pzoop ~                                                               
##     chla              0.550    0.137    4.005    0.000    0.550    0.530
##     corbic           -0.106    0.132   -0.805    0.421   -0.106   -0.102
##     flow             -0.059    0.128   -0.459    0.647   -0.059   -0.059
##     temp              0.149    0.126    1.182    0.237    0.149    0.143
##     secchi            0.091    0.127    0.722    0.471    0.091    0.088
##     hzoop             0.211    0.150    1.401    0.161    0.211    0.197
##   fish ~                                                                
##     hzoop             0.224    0.157    1.433    0.152    0.247    0.237
##     pzoop            -0.259    0.141   -1.837    0.066   -0.285   -0.292
##     corbic            0.113    0.141    0.798    0.425    0.124    0.123
##     flow             -0.021    0.136   -0.153    0.878   -0.023   -0.024
##     temp              0.115    0.136    0.844    0.399    0.126    0.125
##     secchi           -0.353    0.135   -2.611    0.009   -0.389   -0.384
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .estfish           0.000                               0.000    0.000
##    .chla              0.925    0.204    4.528    0.000    0.925    0.948
##    .hzoop             0.587    0.130    4.528    0.000    0.587    0.637
##    .pzoop             0.545    0.120    4.528    0.000    0.545    0.517
##    .fish              0.616    0.136    4.528    0.000    0.746    0.746
## 
## R-Square:
##                    Estimate
##     estfish           1.000
##     chla              0.052
##     hzoop             0.363
##     pzoop             0.483
##     fish              0.254
labelssouth <- createLabels(modfitsouth, cnames)

# residuals(modfitsouth)
# modificationindices(modfitsouth)

Nice plots

Original units

West

North

South

Detrended

West

North

South

Save updated SEM diagrams

plot_grobs <- map(list(figW, figN, figS), ~convert_html_to_grob(.x, 2000))
combined_figure <- ggarrange(plotlist=plot_grobs, labels="auto", 
                             font.label=list(size=9)) %>%
  annotate_figure(top = text_grob("Annual SEM (regional)",
                                  color = "black",
                                  face = "bold",
                                  size = 9))# + 
  #theme(plot.margin = unit(c(1,0,0,0), "lines"))

ggsave('./fig_output/sem_annual_regions.png', combined_figure, width=6, height=6,
       dpi=300)